################################### Packages ##########################################
library(tidyverse)
library(fixest)
library(did)
library(dotwhisker)
library(broom)

################################### Set Working Directory ##########################################

# Getting the path of your current file
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

################################### functions ##########################################
# IHS Transformation
ihs=function(x) log(x+sqrt(1+x^2)) # to revers ihs use sinh
# Demeans and IHS Transforms
idm=function(x) ihs(x)-mean(ihs(x),na.rm=T)
# Demeans
dm=function(x) x-mean(x,na.rm=T)


### main data
up = read.csv("../Data/data_main_final.csv",stringsAsFactors = F)%>%
  ### The following codes replicates the transformations from the main regressions (recreates the regression dataset)
  group_by(stock_id)%>% # for demeaning at the stock level
  arrange(year)%>% # for calculating the lagged variables
  mutate(
    ffmsy = ihs(ffmsy), # The demeaning is only important for the interacted variables, for the non-iteracted variables nothing changes
    gdp=idm(gdp_const),
    gdp2=idm(gdp_const)^2,
    tfp=idm(rtfpna),
    fao_ag_credit = idm(fao_agriculture_credit),
    fao_to_credit = idm(fao_total_credit),
    credit=idm(domestic_credit),
    credit_instrument = idm(credit_instrument),
    interest_factor=dm(1/(1 + lending_interest/100)),
    inflation_factor=dm(1/(1+inflation_cp/100)),
    bbmsy1 = idm(lag(bbmsy,1))
  )%>%
  # The following code is event study specific
  ungroup()%>% # this is very important for the id
  mutate(id = as.numeric(as.factor(paste(stock_id,country,sep="_"))))%>% # creates numeric id
  group_by(id)%>%
  mutate(
    event_time = year - eez_declaration_year, # event time in before and after eez declaration
    max_pr = max(pr,na.rm=T), # maximum treatment level, note that this has to be after the EEZ implementation (also in the main specification the pr level is constant after the EEZ implementation)
    pr_dummy  = ifelse(max_pr>0.5,1,0), # dummy if property rights security exceeds 50% at one point in time, this is for the fixest approach
    eez  = ifelse(max_pr>0.5,eez_declaration_year,0) # this part defines the treatment group: stocks with pr>0.5 are treated and have a treatment year. This is for the did
    )%>% 
  as.data.frame()



################################## Event study: Figure 5: Property rights and resource extraction (continuous treatment) ######################################

### Event study with continuous treatment, baseline in t-1, and some control as main property rights regression
e1 = feols(ffmsy ~ i(event_time, max_pr, ref = c(-1)) + gdp + I(gdp^2) + credit | 
             stock_id + country + year, 
           cluster = ~ country + year + scientific_name,
           up%>%filter(abs(event_time)<=20))


### Event study with binary treatment (one if pr>0.5), baseline in t-1, and some control as main property rights regression
e2 = feols(ffmsy ~ i(event_time, pr_dummy, ref = c(-1)) + gdp + I(gdp^2) + credit|
             stock_id + country + year, 
           cluster = ~ country + year + scientific_name,
           up%>%filter(abs(event_time)<=20)) # 



### creating data frames from regression tables
t1 = tidy(e1)%>%
  mutate(model = "PR continuous")
t2 = tidy(e2)%>%
  mutate(model = "PR binary")

### combine both data frames and clean output for plotting
t = full_join(t1,t2)%>%
  mutate(term = gsub("event_time::|:max_pr|:pr_dummy","",term),
         term = as.numeric(term)
  )%>%
  filter(!is.na(term) & term>=(-20) & term <=20)%>%
  add_row(term = c(-1,-1),
          estimate = c(0,0),
          std.error = c(0,0),
          statistic = c(0,0),
          p.value = c(0,0),
          model = c("PR continuous","PR binary"))%>%
  arrange(desc(term))


### plots output
dwplot(t)+
  theme_bw()+
  coord_flip()+
  geom_vline(xintercept =0,col ="grey50")+
  geom_hline(yintercept = as.factor(0),col = "grey50")+
  scale_colour_manual(values = c( "black","indianred"))+
  theme(legend.position="bottom",
        text=element_text(size=12,  family="serif"))+
  labs(title = "",y = "Event time [Years]",color="",
       x = "Estimate and 95 % confidence interval")+
  scale_y_discrete(breaks = c("-20", "-15", "-10", "-5","0",
                              "5","10","15","20"))

# prints output
ggsave(filename = "../Figures/event_fixest_final.pdf", device = "pdf",width = 20, height =13, units ="cm")




############################# Figure 6: Property rights and resource extraction (binary treatment) using did package by Brantly Callaway and Pedro H.C. Sant'Anna ###################
# number of years
period = unique(up$year)

# creates balanced panel
up_balanced = up%>%
  group_by(id)%>%
  filter(length(id) == length(period))%>% # creates balanced panel 
  select(id, country,year,ffmsy,eez_declaration_year,pr,eez)%>%
  as.data.frame()


## never treated
cs1 = att_gt(yname = "ffmsy",
            tname = "year",
            idname = "id",
            gname = "eez",
            # xformla = ~ gdp + gdp2 + credit, # does not work because of missing gdp and credit observations
            data = up_balanced, 
            control_group = "nevertreated",
            base_period = "varying", # "universal" just increases the confidence intervals of the pre-period - our specification is more conservative
            clustervars = "country",
            allow_unbalanced_panel = F
            
)

summary(cs1)

# event study plot
d1 = aggte(cs1, type = "dynamic",na.rm = TRUE)


# problem:  Not returning pre-test Wald statistic due to singular covariance matrix
# this is the result of the common fisheries policy of the EU that creates colinearity in pr of the same species across eu countries


# average treatment effect
aggte(cs1,type = "simple", clustervars = "country")



#### not yet treated

cs2 = att_gt(yname = "ffmsy",
             tname = "year",
             idname = "id",
             gname = "eez",
             data = up_balanced,
             control_group = "notyettreated",
             clustervars = "country",
             allow_unbalanced_panel = F
             
)

d2 = aggte(cs2, type = "dynamic",na.rm = TRUE)

# average treatment effect
aggte(cs2,type = "simple", clustervars = "country")


#### combine both results in one data frame

t1 = tidy(d1)%>%
  mutate(model = "Never treated")

t2 = tidy(d2)%>%
  mutate(model = "Not yet treated")


t = full_join(t1,t2)%>%
  # full_join(.,t3)%>%
  mutate(term = event.time
         # gsub("ATT","",term),
         # term = gsub("[()]", "", term),
         # term = as.numeric(term)
         )%>%
  filter(abs(event.time)<=20)%>%
  arrange(desc(term))

### Plot results
dwplot(t)+
  theme_bw()+
  coord_flip()+
  geom_vline(xintercept =0,col ="grey50")+
  geom_hline(yintercept = as.factor(0),col = "grey50")+
  scale_colour_manual(values = c("steelblue1","steelblue4"))+ # , "black"
  theme(legend.position="bottom",
        text=element_text(size=12,  family="serif"))+
  labs(title = "",y = "Event time [Years]",color="",
       x = "Estimate and 95 % confidence interval")+
  scale_y_discrete(breaks = c("-20", "-15", "-10", "-5","0",
                              "5","10","15","20"))

### Print results
ggsave(filename = "../Figures/event_did_final.pdf", device = "pdf",width = 20, height =13, units ="cm")





